Which parallel transport for the statistical analysis of longitudinal deformations?

Follow-up imaging studies require the population-wise evaluation of the anatomical changes over time for specific clinical groups. A typical example is the analysis of structural brain changes with aging in Alzheimer’s disease. The longitudinal morphological changes for a specific subject can be evaluated through the non-rigid registration. However, to perform a longitudinal group-wise analysis, the subjectspecific longitudinal trajectories need to be transported in a common reference. The analysis of existing transport methods suggests that the inconsistency of the discrete approximations with composition could drown the small longitudinal deformations into the large inter-subject variability. Instead of approximating a continuous parallel transport in an inconsistent discrete setting, we propose to rely on a carefully designed discrete construction that intrinsically respects all the symmetries on the problem: the Schild’s Ladder procedure. We show that the Schild’s ladder procedure can be implemented extremely efficiently with great stability in the framework of diffeomorphisms parameterized by stationary velocity fields (SVF). The proposed algorithm proves to overcome other methods in terms of consistency with respect to the resampling of most scalar measures of deformations, which opens the way to stable and hopefully more powerful group-wise statistics of longitudinal deformations. 1 Deformations & Morphometry Computational anatomy is an emerging discipline at the interface of geometry, statistics and image analysis which aims at modeling and analyzing the biological shape of tissues and organs. The goal is to estimate representative organ anatomies across diseases, populations, species or ages, to model the organ development across time (growth or aging), to establish their variability, and to correlate this variability information with other functional, genetic or structural information. Following D’Arcy Thompson [4], one often assumes that there is a template shape or image (called an atlas in medical image analysis) which represents the standard anatomy, and that the variability is encoded by deformation of that template towards the shape of each subject or its evolution in time. This approach was pushed forward by Grenander and Miller [7] and turned into a mathematically grounded framework by providing diffeomorphic space deformations with a sufficiently regular right invariant metric [18], leading to the so called Large Deformation Diffeomorphic Metric Mapping (LDDMM) framework. More generally, the local analysis of the deformation retrieved by any image registration algorithm is an increasingly popular method to statistically study differences in brain anatomy. This was coined as Deformation-based morphometry (DBM) by Ashburner [3]. In this context, one should distinguish two types of deformations: the longitudinal analysis of time series of a single subject gives an insight on the development of the subject specific anatomy across time (growth, remodeling, aging), while the cross-section analysis of a group of subjects gives an information on the variability of the organ shape in that population. Distinguishing the two types of deformations is not an easy task and should be performed with care. 1.1 Cross-section designs If the shape variability due to evolution is larger than the intersubject one, then one may consider the inter-subject variability as noise and perform a regression with respect to time even with only one time-point per subject. For instance, the shape of the left and right ventricles of the heart is fairly well defined and only its size normally depends on the age of the child in normal subjects. When we consider pathological case, for instance repaired Tetralogy of Fallot, a remodeling effect perturbs the normal growth and drastically affect the shape of the right ventricle. In that case, a partial-least-squares regression of the shape changes with respect to age allows us to model the principle effects of the remodeling of the right ventricle [12]. The resulting statistical model of growth turns out to have an anatomically meaningful interpretation which could be assessed by correlating its shape with other clinically important variable. For instance, the severity of regurgitation is associated with a bulging of the outflow tract, a dilatation of the apex and a global RV dilatation. Such a linear (geodesic) deformation model in the LDDMM framework can be extended to a geodesic by part regression when more complex deformation patterns are sought, potentially including a time-warp independent of the spatial deformation. When applied to different species (bonobos vs chimpanzees) or to diseases (autism vs control), this model suggests that the change in the speed of evolution might sometimes be more important than the shape differences [5]. 1.2 Longitudinal designs In most cases, longitudinal changes of organs are quite small compared to inter-subject variations. In pathologies such as dementia and in particular Alzheimer’s Disease (AD), it is quite clear now that longitudinal measures are more sensitive (but way smaller) than cross-sectional differences for distinguishing normal aging from mild-cognitive impairments (MCI) from AD. As longitudinal measures are done at the individual level, we need to normalized them in a standard reference anatomy in order to obtain statistical significant localization of these changes at the population level (group-wise statistics). In many cases, one consider a scalar summary of the longitudinal changes fi(x, t) at the subject-specific level and a (diffeomorphic) deformation Ψi(x) relating the template to the inter-subject coordinate systems. In fMRI, the scalar function can be the z-score of the statistical correlation of the input to the signal. For morphological changes, the map of the determinant of the Jacobian fi(x, t) = |∇φi(x, t)| of the longitudinal deformation φi(x, t) represents the local volume change. Other measurements are possible: the log of the above determinant can be interpreted as the local flux of the deformation through the volume element [9], the L2 norm of the displacement field, the local Henky or St Venant-Kirchoff elasticity energy [13], the norm of the momentum in LDDMM methods, etc. In all these cases, the natural normalization method is to resample the function values in the template coordinate space f̂i(x, t) = fi(Ψi(x), t). Notice however that when the map considered is a density (e.g. gray or white matter density in VBM), a volume-preserving modulation is applied to correct for the volume change of the integration element [2]. This is the simplest example of the more complex situation where the geometric nature of the longitudinal measure should be taken into account. For instance vector fields should be reoriented according to the gradient of the inter-subject deformation in order to take into account the local rotation and skew of the deformation. Likewise, it is considered in diffusion tensor imaging that the directions of diffusion should be reoriented, but not the eigenvalues which have a physical meaning. This lead to potentially multiple ways to define the optimal rotation of a tensor, like Preservation of the Principal Directions (PPD) or the Finite Strain (FS) [16]. 2 Normalizing longitudinal deformations Transporting a scalar summary of the changes over time is numerically stable and allows elaborating statistical difference maps at the group level. However, this does not allow to go much beyond this detection level. In particular, one cannot model the ”average” group-wise deformation and its variability, nor transport it back at the subject level to predict what will be the future deformation. To reach such a generative model of the longitudinal deformations, we should normalize the deformations themselves and not just some of its components. This involves defining a method of transport of the longitudinal deformation parameters along the inter-subject change of coordinate system. It is important to recall here that the nature of the longitudinal and inter-subject deformations could be very different, Thus, in terms of geometry, there is a-priori no reason to consider the same geometric structure (connection or metric) for both deformation spaces. Depending on the considered parameterization of the transformation (displacement fields, stationary velocity fields, initial momentum field...), different approaches have been proposed in the literature. A simple method consists in considering the longitudinal transformation parameters as vectors fields which are reoriented by the Jacobian matrix of the inter-subject mapping. A clear drawback of this method is that it is only valid for small longitudinal displacements and no consistency can be expected along the temporal trajectory. A more elaborated method consist in assuming that the longitudinal deformation occurs in the subject specific coordinate system, and that the inter-subject deformation is a (static) change of coordinate system which is valid for all time frames. In that case, the trajectory in the template space will be defined by the conjugate action φ̂i(x, t) = Ψ(φi(Ψi(x), t)). This transformation conjugation was used for instance in [15, 5]. In the LDDMM framework, [17, 14] proposed to to parallel transport the initial momentum of the longitudinal deformation along the inter-subject geodesic diffeomorphism. The proposed method is fully consistent with the chosen right-invariant metric on diffeomorphisms, but unfortunately requires that the same metric is used for both longitudinal and inter-subject deformation. In almost all cases, one can question the numerical stability of the scheme and the consistency of the discrete encoding of transformations with the composition of transformations. In particular, the large inter-subject deformations may easily induce large approximation errors on the transport of the small longitudinal deformations. 2.1 Parallel transport using Schild’s ladder Instead of defining properly the parallel transport in the continuous setting and approximating it in an inconsistent discrete setting, we proposed in [11] to rely on a carefully designed discrete construction that intrinsically respects all the symmetries on the problem: the Schild’s Ladder procedure algorithm, initially introduced in the last century in the field of the general relativity by the physicist Alfred Schild [6]. Maybe even more importantly, this procedure allows to parallel transport vectors along any curve and not just geodesics. The method is based on the construction of a geodesic parallelogram illustrated in Fig.1. To transport vector A at point P0 to point P1 along the curve elementC of length , we first build a point P ′ 0 = expP0( A) in the direction of A. Then we compute the midpoint P2 = expP ′ 0( 1 2 logP ′ 0(P1)) on the geodesic from P ′ 0 to P1. The geodesic starting from P0 whose midpoint is P2 is the second diagonal of the geodesic parallelogram and defines the point P ′ 1 = exp(2 logP0(P2)). The parallel transport of vector A at point P1 is A′ = logP1(P ′ 1). The whole procedure is then iterated along the successive elements of the curve to parallel transport a point along the whole curve. FIG. 1: Right: The geodesic parallelogram forming each step of the Schild’s Ladder. Left: the procedure is iterated along the curve to constitute the steps of the ladder. The Schild’s ladder procedure only requires the computation of (Riemannian) exponential and logarithms, and thus can easily be implemented for any manifold provided that we have these basic algorithmic bricks. In fact, one can show that the procedure is valid for any affine connection space provided that the connection is symmetric [8]. 2.2 Implementation for SVFs acting on images The Schild’s ladder procedure can be implemented particularly simply in the framework of stationary velocity fields (SVF), which generate one-parameter subgroups of transformations. This infinite-dimensional analogous of the Log-Euclidean framework for tensors and linear transformations was introduced in [1] to parameterize diffeomorphic deformations as the flow a stationary velocity field φ = Exp(u) (instead of a time-varying one in LDDMM). A key point is that the scaling and squaring algorithm allows computing the exponential in a logarithmic time with respect to the number of time steps normally needed to integrate the ODE. As the inverse deformation is parameterized by the opposite of the SVF, we can very easily enforce inverse consistency. Moreover, the Baker-Campbell-Hausdorff (BCH) formula gives a very efficient way to approximate the composition of two deformations directly in the log-domain. From the theoretical point of view, the SVF framework can be justified by dropping the Riemannian metric for the more general setting of affine connections. Indeed, there is a unique symmetric connection which is both left and right invariant on a Lie group: the canonical symmetric Cartan connection. The geodesics of this connection are the left (and right) translates of the one parameter subgroups. Thus, SVFs are parameterizing the geodesics starting from identity. Notice that the space is not flat with this connection as we have R(X,Y )Z = −1/4[[X,Y ], Z] for any left invariant vector fields X,Y, Z. Of course, the bi-invariance (which also implies the equivariance by inversion) cannot be obtained for free. We loose in general the property that any two points can be joined by a smooth geodesic, which was automatically verified in Lie groups with right invariant metrics (e.g. LDDMM). In infinite dimension, the exponential maps might even not be locally diffeomorphic. These points remain to be carefully investigated but in practice the set of diffeomorphisms that are reachable by one parameter subgroups seems to be sufficient large to encompass all reasonable anatomical deformations. Last but not least, we are interested here into deformation of images. Thus, we should replace the points in our space by images, and consider that the tangent vectors are the SVFs that drive the deformation of the current image: expI(v) ' exp(u)?I . Likewise, the log of image J with respect to image I is the SVF that best registers I to J . In this context, the Schild’s ladder step transporting the SFV u at image I0 to the (template) image T0 can be rewritten as follows (Fig. 2). First, we warp Image I0 according to the SVF u to get I1 = Exp(u) ? I0. Second, we register image I1 to image T0. The result is the SFV v and the midpoint image would be I1/2 = Exp(v/2) ? I1 or symmetrically I1/2 = Exp(−v/2) ? T0. Finally, we get by